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' We derive exact solitonic solutions of a class of Gross-Pitaevskii equations with time-dependent 
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harmonic trapping potential and interatomic interaction. We find families of exact single-solitonic, 
multi-solitonic, and solitary wave solutions. We show that, with the special case of an oscillating 
trapping potential and interatomic interaction, a soliton can be localized indefinitely at an arbitrary 
I position. The localization is shown to be experimentally possible for sufficiently long time even with 

only an oscillating trapping potential and a constant interatomic interaction. 
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^ I. INTRODUCTION 

The experimental realization of solitons in Bose-Einstein condensates [ll, S 0j S @] has stimulated intense interest 
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in their properties [aJSUa llfl, HH, |1J, OJ, OJ. 
the old [TgI. IITI. [l8l.ll9n20| and more recent 



15 1 . The inhomogeneity provided by the trapping potential has renewed 
21| interest in the different aspects of one- and multi-solitons dynamics 



^ in inhomogeneous potentials. 

In the experiments of Refs. [1, stable bright solitons were created and set in a particle- like center-of-mass motion. 
I t , The wave nature of solitons was revealed when two adjacent solitons repelled each other as a result of their phase 
Ci ' difference On the other hand, it is established that bright solitons collapse when the number of atoms exceeds a 
H . certain limit (2^ . This results from the attractive Hartree energy overcoming the repulsive kinetic energy pressure. 



I • One of the methods proposed to stabilize the soliton against collapsing is to rapidly oscillate the interatomic interaction 



or the trapping potential [23[. Obviously, the inhomogeneity imposed by the trapping potential plays an important 
role on the dynamics and stability of solitons. 
^ ' The evolution of solitons is approximately described by the inhomogeneous nonlinear Schrodinger equation known 
I— as the Gross-Pitaeviskii equation p4, 25]. The approximation stems from the fact that the Gross-Pitaevskii equation 
is a mean-field approximation of the exact iV-particle Schrodinger equation. While in some cases the soliton dynamics 
<^ , obtained by these two equations disagree the Gross-Pitaevskii equation often gives accurate results. Theoretical 

■ studies performed to account for the observed behavior of solitons were conducted by solving the Gross-Pitaevskii 
[ equation with numerical, perturbative, or variational methods. Much less effort was devoted to finding exact solu- 

(SJ ■ tions of this equation [ll, [l^, [H, [H, [13, [Ml, [H, [H, [s^l- In addition to providing rigorous insight, exact solutions 
acquire valuable importance when problems such as soliton-soliton collisions and soliton interaction with potentials 

ly-j are addressed. In such cases, formulae for the force between solitons or soliton's effective mass can be derived 'ss']. In 
addition, obtaining such exact solutions allows for testing the validity of the Gross-Pitaevskii equation at high soliton 

On densities and obtaining the long-time evolution of the soliton where numerical techniques aught to break down. 

■ Here, we further explore exact solitonic solutions of the Gross-Pitaevskii equation. Specifically, the goal of this paper 
^ I is two-fold. First, we investigate the existence and properties of solitonic solutions in the presence of time-dependent 

trapping potential and interatomic interaction. Secondly, we focus on the effect of an oscillating trapping potential 
and interatomic interaction on the center-of-mass motion of the soliton. We consider here only the case of attractive 
^ . interatomic interactions which allows for bright solitons. 
■ ■ ' The first goal is achieved by employing the Darboux transformation method ^36^ to derive families of exact solitonic 
solutions of a class of Gross-Pitaevskii equations. It should be noted that the main solution we derive here (Eq. ([211)), 
which corresponds to a harmonic expulsive potential, reproduces a special case of the general solution found by Serkin 
et al. corresponding to a combination of a harmonic and linear potentials [2l)]. Hence, the significance of the first 
goal is mainly presenting a systematic method of generating exact solutions. 

For the second goal, we considered the special case of oscillating strengths of the harmonic trapping potential and 
interatomic interaction. Interestingly enough, it turns out that such oscillations not only stabilize the soliton against 
shrinking, but also make it possible to localize it at an arbitrary position. This is the main result of this paper. The 
possibility of localizing the soliton is then discussed from an experimental point of view. To that end, we considered the 
nonintegrable, though experimentally simpler, case of only an oscillating trapping potential and constant interatomic 
interaction. Here too, the soliton can be localized, though not indefinitely as before. The long-time evolution shows 
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that the sohton continues to be trapped but will be oscillating around the minimum of the harmonic potential. For a 
typical experimental setup, we show that the soliton can be localized around its initial position for a time period of the 
order of, or even larger than, the lifetime of the soliton. The soliton localization suggests a management mechanism 
for the soliton position and speed that may have applications in various situations such as soliton-soliton collisions 
and soliton interaction with potentials. 

The rest of the paper is organized as follows. In the next section, we present the general form of the Gross-Pitaevskii 
equation to be solved. In section Hill we use the Darboux transformation method to derive the new solitonic solutions. 
We then discuss their properties, dynamics, and localization. In section HVl we discuss the experimental feasibility of 
realizing soliton localization. We end in section IVl with a summary of our main results and conclusions. 



II. THE GROSS-PITAEVSKII EQUATION 



The Gross-Pitaevskii equation describing a Bose-Einstein condensate trapped by an axially-symmetric harmonic 
potential with attractive interatomic interactions is given by 
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where is the absolute value of the s-wave scattering length, and uj^ and uj± are the characteristic frequencies of the 
harmonic trapping potential in the axial and radial directions, respectively. 

When the confinement of the Bose-Einstein condensate is much larger in the y and z directions compared to 
the confinement in the x direction, the system can be considered effectively one-dimensional along the x direction. 
The three-dimensional Gross-Pitaevskii equation can then be integrated over the y and z directions to result in a 
one-dimensional Gross-Pitaevskii equation 0, [s^l 



^rl—i}{x,t) = 
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where A — uj^juj^. Scaling length to = ^ %/ mcux , time to 1/uJx, and il}{x,t) to l/y/2Xax, the Gross-Pitaevskii 
equation takes the dimensionless form 



'\§^ + \p{t)x^ -aq{t)\^{x,tf 



ip{x,t), 



(3) 



where a = as/a^ is the scaled scattering length. The dimensionless general functions p{t) and q{t) are introduced to 
account for the time-dependencies of the strengths of the trapping potential and the interatomic interaction. 

In the subsequent section, it is shown that this equation is integrable only if p{t) and q(t) are parametrically related 
as follows: p{t) = 7(t) — 7(i)^ and q{t) = exp(7(i)), where ^{t) is an arbitrary real function (see Eq. ([7])). We find 
exact solitonic solutions of Eq. ^ below in terms of the function "f{t). A host of rich and interesting physical systems 
are described by such a class of solutions. In particular, an oscillating form of j{t) will be considered in this paper. 



III. DARBOUX TRANSFORMATION AND THE EXACT SOLUTIONS 



A. Darboux Transformation 



Applying the Darboux transformation method on nonlinear partial differential equations requires finding a linear 
system of equations for an auxiliary field ^{x,t). The linear system is usually written in a compact form in terms of a 
pair of matrices as follows: —TJ-^ and = The matrices U and V, known as the Lax pair, are functionals 

of the solution of the differential equation. The consistency condition of the linear system "^xt = "^tx is required to be 
equivalent to the partial differential equation under consideration. Applying the Darboux transformation, as defined 
below, on ^ transforms it into another field For the transformed field '4'[1] to be a solution of the linear system, 

the Lax pair must also be transformed in a certain manner. The transformed Lax pair will be a functional of a new 
solution of the same differential equation. 

Practically, this is performed as follows. First, we find the Lax pair and an exact solution of the differential equation, 
known as the seed solution. Fortunately, the trivial solution can be used as a seed, leading to nontrivial solutions. 
Using the Lax pair and the seed solution, the linear system is then solved and the components of '4' are found. The 
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new solution is expressed in terms of these components and the seed solution. The following detailed derivation of 
the new solution clarifies this procedure further. 

Using our Lax Pair search method [11], we find the following linear system which corresponds to the class of 
Gross-Pitaevskii equations we are interested in: 

*^ = CJ * A + U (4) 
*4 = • * • A • A + C («U + a; -/{t) J) •* ■ A + V ■ *, (5) 

where, 

VPfr f)-fM^,t) M^,t)\ f 1 

Ui(a;,i) M^,t) r I -1 



V 



Ai \ / V^Q(x,t) 

A2 J' 1^ -^^g*(x,^) 

ia\Q{x, i)| 72 ^\x^Q{x, t) + i^Q.^{x, t)/2 

-y^XxjQ*{x, t) + i^Q%{x, t)l2 ~ia\Q{x, t)\^ jl 



— exp(7(i)), and Ai and A2 are arbitrary constants. For convenience, we presented the matrices in 
terms of the function Q(x,t) which is related to the wave function as follows Q{x,t) — t)e'^'*''^*-'+*^'^*^^ 

It should be emphasized that while applying the Darboux transformation is almost straightforward, finding a linear 
system that corresponds to the differential equation at hand is certainly not a trivial matter. Usually, this is found by 
trial and error, or by starting from a certain linear system and then finding the differential equation it corresponds to. 
In Ref. [385 , we have introduced a systematic approach to find the linear system which we describe here briefly. The 
partial derivatives of the auxiliary field, and ^t, are expanded in powers of A with unknown matrix coefficients. 
The expansions are terminated at the first order for and the second order for since this will be sufhcient to 
generate the class of Gross-Pitaevskii equations under consideration. The higher order matrix coefficients turn out to 
be essentially determined by the zeroth order matrix coefficients U and V. To find the matrices U and V, we expand 
them in powers of the wavefunction -0(2;, i), its complex conjugate, and their partial derivatives. The coefficients of 
the expansions are unknown functions of x and t. Substituting these expansions in the consistency condition (Eq. ^ 
below) we find a set of equations for the unknown function coefficients. Finally, by solving these equations the Lax 
pair and consequently the linear system will be determined. The linear system found here is a generalization to that 
of Zakharov-Shabat for homogeneous Gross-Pitaevskii equation [40|. 

For to be a solution of both Eqs. (|4]) and (O, the consistency condition "i^xt = ^tx must be satisfied. This 
condition leads to the following relation between the matrices U and V 

Ut- V, + [U,V] =0, (6) 

where [U, V] is the commutator of U and V. Substituting the above expressions for U and V in the last equation, 
we obtain the Gross-Pitaevskii equation 



dt 



0(x,t), (7) 



and its complex conjugate. This equation shows that the functions p{t) and q{t) are parametrically related to each 
other through the general function 7(t). In Ref. 21], Serkin et al. solve a nonautonomous Gross-Pitaevskii equation 
that is similar to Eq. ([7]) but with additional time-dependent dispersion and linear potential. Similar to our conclusion 
[38t , the authors of this reference and Ref. [13] found that relations between the coefficients must be obeyed for the 
the model to be integrable. We focus here on the specific case of only a harmonic trapping potential and constant 
dispersion. The linear system of 8 equations, Eqs. (|3|) and read explicitly 

V'lx - </'iVieM»^-'+^)V;o - V2Ai0ie^ = 0, (8) 
- 02\/^eM^i'-'+^)v;o - V2A202e^ = 0, (9) 
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(^1, + V2Xi e>i + 7A*V'iv^eKT-"'^) = 0, (10) 

hco + V2X2 e>2 + Vo*V'2A/^eK^-"'^) = 0, (11) 

V'lt - «V'ie^ (^2eTA? - i\/2a;7Ai + ^ [V-o p a) - ^(^1 V^eK^^^'+'z) (^^^.^^ + (2^V2e^Al + x-r)) = 0, (12) 

V'2t - ^^^26^ (^2e^A^ - iV2a;7A2 + ^ [V-o |' ~ i</'2V^eH*^"'+'') (i^o. + (2iV2e^X2 + x-r)) = 0, (13) 

01* + ^V'lV^eK^-"'^) (v^* (2iV2e'^Ai + 2:7) - ^V^.) + i(/)ie^ (i jV'o |' a + 2Ai (2ie^Ai + 722:7)) = 0, (14) 

02* + ^?/'2\/^eK^-"'^) (2iV2eTA2 + X7) - i^S.) + ^026^ (i |' a + 2A2 (2ie''A2 + 722:7)) = 0, (15) 

where ipQ^Xjt) is an exact seed solution of Eq. ([7|). These equations reduce to an equivalent system of 4 equations 
with nontrivial solutions by making the following substitutions: Ai — — A2, '02 = 0i, and 02 = —ipi- Using the trivial 
solution, ^pQ(x,t) = 0, as a seed, the linear system will have the solution 

1pl(X,t) = cie , (16) 

(p2(X,t) ^ 026 , (17) 

where ci and C2 are real arbitrary constants of integration. 

Consider the following version of the Darboux transformation [36[ 

*[1] = A-cr*, (18) 

where '4'[1] is the transformed field, a — 'i'o • A • 4*0"^ ■ Here 'i'o is a known (seed) solution of the linear system 
Eqs. ((8]fT5|. For the transformed field ^[1] to be a solution of the linear system, the matrix U for instance must be 
transformed as [s^ 

U[l] = (7 • U • (7~^ + CT^ • (T-\ (19) 

where is the inverse of cr. This equation gives the new solution in terms of the seed solutions of the Gross-Pitaevskii 
equation, ip(){x,t), and the linear system, 4'o, which reads 

^(x,t) = ^o(x,t) + ^(Ai + Xl)e-^^W^'/^+'rW/2^^t^ (20) 

B. Exact single-solitonic solution 

Using 4'o{x,t) = and Eqs. and ^T7} . the new solution takes the form 

4V2Ai^ciC2^ e^W/2+.^(*) 

where 9{t) = -4 (Ai,^ - Ai^^) / e^-i-W dt+x (2 e-^W Ai» - 2; 7(i)/\/8) , /3(t) = 2V2Xire'^^*\ and 2;o(t) = 
8 Ah Air / e^-^Wdt. Here ci, C2, Ai, and A2 are arbitrary constants. The subscripts r and i denote real and imaginary 
parts, respectively. Substituting ci = exp((^i), C2 = exp((52), where Si and (^2 are arbitrary constants, completing the 
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square in the phase factor, and normahzing ip{z,t) to this solution can be recast in the following more appealing 
form 



where 



= ^^^X^e'^(*)/2+^^(-*)sech(^^e-W(x-Xe„.(i)) 



(22) 

(23) 
(24) 



C3 



it? + 7^{-MNf + («o + 2:0 7(0))' g{t) 



(25) 



f;(i) = e^'*''^* ^dt'. The constant C3 corresponds to an arbitrary overall phase. 

This solution corresponds to a sech-shaped soliton containing N atoms with a time-dependent center-of-mass Xcm(^) 
and time-dependent width 2 exp(— 7 1) /N a. The linear part of the phase profile shows that the soliton is moving with 
a center-of-mass velocity v{t) — icm(i)- The quadratic part corresponds to a phase chirp associated with the quadratic 
trapping potential. 

The simple choice ^{t) = constant^ which corresponds to the homogeneous Gross-Pitaevskii equation, gives Xcra{t) = 
vqI + xq. The solitonic solution in this case corresponds to the well-known sech-shaped soliton with a center-of-mass 
moving with a constant velocity vq and starting the motion at the initial position xq. In the limit ^(t) — *■ 0, the 
solution, Eq. (j^ . reduces to 



tp{x,t) = 



(*)sech 



Na 



[x - (xo + Wot)) 



(26) 



where (j){t) = C3 -I- ((aiV)^ -I- 4uq) t/8 + V[){x — (xq -f fo 0)- This is the well-known sech solution of the homogeneous 
Gross-Pitaevskii equation. For 7(t) = constant x t we get a Gross-Pitaevskii equation with an expulsive harmonic 
potential and exponentially growing interatomic interaction [30]. In principle, we can choose any form of 7(i), but one 
should keep in mind that the interatomic interaction strength is proportional to exp(7(i:)). Such a time-dependence 
may not be realistic from an experimental point of view for any 7(t). 

In Ref. [m Serkin et al. have already derived the exact solitonic solution of a Gross-Pitaevskii equation with a 
linear and quadratic potentials and time-dependent dispersion 
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l-g^V{X,t) = 



Djt) n^t) 



x^ - 2a{t)x - a R{t)\il){x, t)[' 



■ip{x,t). 



(27) 



The authors show that this equation is integrable o nly if the functions D{t), ri(t) and R{t) are related to each other 



through the integrability condition (Eq. (2) in Ref. 21[). The special case of only a quadratic potential and constant 
dispersion corresponds to the Gross-Pitaevskii equation considered in this paper. Substituting a{t) = and D(t) — 1 



in Eq. (2) of Ref. [21], the integrability condition simplifies to fl 



7, which shows that, in this special case, 



the previous equation is indeed equivalent to Eq. ([7]). Therefore, substituting a{t) = and D{t) = 1 in the general 
solution of Eq. (|27p. namely Eq. (8) in Ref. \21\, should result in our solution, Eq. (|^^ . It turns out, however, that 
the two solutions do not match exactly. The solution of Ref. [2l[ corresponds to a soliton located initially at x = 
while in our case the soliton is located initially at the arbitrary position xq . This can be clearly seen by substituting, 
without loss of generality, r]{t) = N a exp(7(t))/4 and K{t) — —vq exp(— 7(0))/2 in the argument of the sech function 
of Eq. (8) of Ref. 0. This results in a center of mass coordinate Xcm{t) = v^git) e~"^^^'>~"^^^\ Since g{Q) = 0, this 
shows that Xcm(O) = in contrast with our case a;cm(0) = Xq. 

The dynamics of the soliton is readily given by Eq. ([M]) . An equation of motion for the center-of-mass Xcm{t) can 



be derived from the Euler's equation of the lagrangian L[x, 
functional is given by 



cm 1 ^cmj 



J^^idxij*dij/dt - E[xcm,Xcm]- The energy 



E[x, 



cm 1 -^cmj 



rix,t) 



29x 



1 



+ 0(7(0-7(0')^'- o«e^'*'lV'(2;,i)l2 



^{x,t), 



(28) 



which results in the equation of motion 



,(0 + (7-7')a;cm(i) =0. 



(29) 

This equation shows that the center-of-mass motion is determined by the function ^{t). Thus, interatomic interactions 
do not affect the center-of-mass motion which is a manifestation of Kohn's theorem |4l| . 
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C. Oscillating trapping potential 



In this section, we consider an oscillating form of "f{t) which results in a trapping potential and interatomic interac- 
tion with oscillating strengths. A first simple choice for j{t) would be for instance cos {ojt + 6). However, in this case, 
the interatomic interaction, which will be proportional to exp (cos ujt + S), oscillates nonlinearly. Such time-dependent 
interatomic interaction may not be possible to realize experimentally. Instead, we use the form 

7(i) = ^ (ai + "2 cos {ujt + S)), (30) 

where ai, a2, lj, and S are arbitrary dimcnsionlcss constants. In this case, the Gross-Pitaevskii equation, Eq. ([7]), 
takes the form 



dtp(x, t) 



2 4 



a2Uj'^ {cos {cut + 6) + sin {cut + - a e("i+"^ |2 



i^{x,t). (31) 



The advantage of this particular form of 7(i) is that, for a2 ^ 1, the amplitude of the oscillation in the interatomic 
interaction can be set to an arbitrarily small value such that the strength of the interatomic interaction can be 
considered practically as constant. Substituting this expression for ^{t) in Eq. (|^5)) . we get 



-a2 ^2 COS {ujt + 5) + a2 sin {ijjt + 5)^^ Xcm — 0. (32) 
The general solution of this equation is readily given by Eq. p4|l . which now takes the form 

The first term of this equation corresponds to a bounded oscillation, but the presence of the integral in the second 
term makes Xcra{t) unbounded. Due to this term, the soliton will be expelled out of the trap, i.e., Xc-m{t oo) = ±oo. 
Therefore, the soliton can be localized, by choosing the parameters such that the prefactor of the unbounded term 
vanishes, namely 

Wo — 2 0^2 sin^. (34) 
With this condition, the center-of-mass of the soliton is given by 

Xc^{t) = a:oe^(™"*-™^ (35) 

This is one of the main conclusions of this paper. It shows that the soliton can be localized at an arbitrary position 
by oscillating the strength of the trapping potential and the interatomic interaction. Without such oscillations, the 
soliton will be expelled out of the trap. In fact, this result also holds for any bounded type of oscillations as can be 
inferred from Eq. Taking vq + Xo7(0) = 0, this equation gives Xcm{t) = xq e'''^^^"'''^'^ which is bounded for any 

bounded 7(t). The different cases of soliton localization and delocalization, described by Eq. ((33)) . are shown in Figs.[T] 
andm In Fig. [1] the soliton is shown to be expelled out of the left (right) side of the trap when vq < xoa2UJ sin (S) /2 
(vq > xoa2U) sin((5)/2), while for vq = xq a2 o; sin(5)/2, the soliton remains localized around its initial position xq. 
In Fig. [21 this is shown with the trajectory of the center-of-mass of the soliton. 

Notice that in order to localize the soliton, no condition on lu was required. Therefore, one may argue that by 
taking uj arbitrarily small, we get a localized soliton in an almost stationary expulsive harmonic trap. This of course 
contradicts the fact that in an expulsive harmonic trap, solitons are expelled away from the center. However, one 
should keep in mind that with our special form of 7(t), namely Eq. (j30[) . the strength of the trapping potential will 
be proportional to a;^. Therefore, a very small value of u corresponds to a shallow potential that approaches the 
homogeneous case for uj — 0. The fact that the strength of the harmonic potential depends on oj leads to conclude 
that, for deep traps, larger trap oscillations frequency are needed to localize the soliton contrary to the case with 
shallow traps where the soliton can be localized with smaller trap oscillations frequency. 

Taking the time average of the trapping potential (w/27r) J^^^'^ dt {■y{t) — j{t)'^)x'^/2 — —a^bj^x^jX^ shows that 
the soliton spends on the average more time in the expulsive trapping potential. One may thus conclude that after 
sufficiently long time the soliton will be expelled out of the trap which contradicts our previous localization result. 
A careful examination of the dynamics shows that this conclusion is incorrect. In spite of the fact that the soliton 
spends more time in the expulsive trap, the inhomogeneity of the trapping potential can compensate for the time 
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difference. The strength of the trapping potential has in general two periods ti and r2, as shown in Fig. [S] The 
period ti depends on 02 and T2 — 27r/w. The strength of the trapping potential is positive for a period of Ti and 
negative for T2 — ti > ti for all 012. Assuming the soliton started the motion at > from rest, it will experience 
at first a harmonic trapping potential for time ti and therefore will move to the left (region of lower potential) for a 
certain distance. At time ti, the trapping potential becomes expulsive for time period T2 — ti and the soliton starts 
to return back. However, it starts now from a point of lower potential than at xq which means that it will experience 
a weaker trapping potential. Therefore, in order to reach the starting point, xq, it needs more time compared to the 
forward part of the motion. If that time matches T2 — ti, the soliton returns back to xq with zero velocity and the 
cycle repeats leading to trapping of the soliton, which corresponds to the middle curve of Fig. [21 On the other hand, 
if the soliton reaches a point a; > a;o , it will be eventually expelled out of the right side of the trap corresponding to 
the upper curve of Fig. [21 and if it reaches x < Xq, it will be expelled out of the left side of the trap corresponding to 
the lower curve of Fig. [51 

To further understand this trapping mechanism, we calculated the trapping potential felt by the soliton through its 
trajectory, namely V{xcm{t)), which is plotted in Fig.[4l In this figure, the dynamics of the soliton is represented by a 
point moving on the curve with a direction that is indicated on each curve. The center-of-mass motion of Fig. [2lcan be 
extracted by tracing Xcm{t) while the point moves along the potential curves. In Fig[3^, the soliton starts at xq = 50 
with initial velocity vq = 1.2 012X01^ sin((5)/2. In this case, the soliton gets drifted by time towards larger values of 
Xcm corresponding to the upper curve of Fig. [1] In Fig. [31;, the initial velocity is vq = a2 xqlu sin (S) /2, which results 
in soliton localization. In this subfigure, the soliton oscillates between Xcm = 50 and Xcm = 61 corresponding to the 
middle curve of Fig.[21 In Fig.[l^, the initial velocity is vq = 0.8 a2 xq lj sin (6) /2 leading to a drift towards lower values 
of corresponding to the lower curve of Fig. [21 To make an analogy with a classical particle moving in a time- 
independent potential, we defined VesiXc-mit)) = V{Xc-m{t))-[V{Xcrni'^/i^))-V{XcniiO))][Xcm{t)-Xcm{0)]/[Xcmi''T/uj)- 

2;cm(0)]. In Figs. [3}d, d, f, we plot T4ff(a;cm) that corresponds to Figs. [1^, c, e, respectively. The dynamics is now 
simplified to that of a classical particle oscillating in a ladder of parabolic potentials. For the nonlocalized soliton, the 
potential minimum is shifted by time to the right (Fig. [41d) or to the left (Fig.[4f). For the localized soliton case, the 
minimum of the potential is stationary. 

D. Another family of exact solutions 

Using the exact solution found above as a seed solution, the Darboux transformation generates a two-solitons 
solution. This kind of solution is useful for studying soliton-soliton interactions which is left for future work. Another 
family of more complicated exact solitonic solutions can be obtained by using a nontrivial seed solution as shown in 
Appendix A. Substituting for j{t) in Eq. (|A2[) . we get the second class of exact solutions. This family of solutions is 
more complicated than the above single-solitonic solution since it involves, in addition to the single-solitonic solutions, 
multi-solitonic and solitary wave solutions. Here, we present briefly the main properties of the solutions of this kind. 
In Fig. [51 we plot the density of a single-solitonic solution showing that the soliton is being expelled out of the center of 
the trap. The oscillation in the trajectory of the soliton is due to the oscillating trapping potential. The discontinuous 
appearing of the soliton is due to the interaction with the background. This trajectory can also be extracted from 
the general solution, Eq. (|A2|1 . by considering the term (Pcf e^"''"^^'''^ -I- rc|e~^"''"'"^^' '^) in the denominator. At 
the soliton's density peak this is the dominant term that determines the position of the soliton. Specifically, the 
position of the peak is given by the condition 2ari — y^A^a; = 0. Using this condition to plot x versus t in Fig. [6l 
we obtain a curve that is identical to the soliton trajectory in Fig. [5l The mean slope of this curve is proportional to 
a/2Ar. Hence, the rate at which the soliton leaves the center of the trapping potential can be delayed by choosing 
the parameters and the arbitrary constants such that a/2Ar is small. For the special case of a = the center-of-mass 
of the soliton will be localized at a; = indefinitely. In this case, the oscillating trapping potential results only in 
oscillations in the width and peak density of the soliton. This is also shown in Fig. [3 where we see a multi-solitonic 
solution with central soliton being localized at a; = and off-central ones oscillating around their initial positions. In 
Fig. [51 we show that for some values of the parameters the dynamics of the peak soliton density can be so drastic 
such that the soliton disappears in the background and reappears at regular discrete times. 

IV. NUMERICAL SOLUTION AND EXPERIMENTAL REALIZATION 

As we have seen in section IIII C| there is a possibility to localize the soliton by oscillating the trapping potential 
and the interatomic interaction in the manner described by Eq. ([?T|l . Such synchronized oscillations may not be 
possible to realize experimentally. Instead, a setup with an oscillating strength of the trapping potential and constant 
interatomic interaction may be experimentally more favorable. This situation can be obtained in our case with the 
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condition Q!2 ^ 1 resulting in the following Gross-Pitaevskii equation 

92 1 



dii){x,t) 



dx^ - cos {cot + 5)x^ -a \ifj{x, t)f 



^P{x,t) (36) 



where J7 = w y/a2/2 and we have set ai — 0. We solve this equation numerically using the exact solution of the 
homogeneous case, namely the solution of Eq. p6p with f2 = 0, as the initial wavefunction. This can be obtained from 
Eq. simply by substituting ai = a2 = 0. The soliton's center-of-mass trajectory is extracted from the resulting 
numerical solution and then plotted versus time as shown in Fig. [5) The trajectory is shown with the filled circles for 
io = 2 and fl — 0.44 and with the empty circles for w = 2 and = 0.14. It is clear from this figure that with smaller 
amplitude of the oscillating trapping potential, the soliton will be locafized for longer periods. The sofid and dashed 
curves show the corresponding trajectories in the presence of the oscillating interatomic interaction as described by 
Eq. ([55]) . The difference between the solid and filled circles curves shows the important role played by the oscillations 
in the interatomic interactions in stabilizing the soliton. On the other hand, the overlap between the dashed and open 
circles curves shows that the effect of interatomic interactions is minor for smaller amplitudes of the oscillation in the 
trapping potential. 

The exact center-of-mass dynamics of the solitonic solution of Eq. (|36p is dictated, according to Kohn's theorem 
[4l| . by the potential — 5^^^ cos {uit + S) x^ independently from the interatomic interaction. Taking advantage of this 
fact, the equation of motion of the center of mass follows 

Xcm{t) ~ fl'^ cos {ujt + S) Xcm{t) ^ 0- (37) 

The general solution of this equation is a linear combination of the sine and cosine Mathieu functions 

Xcmit) = ci C(0, a2, {lj t + 6)/2) + 5(0, aa, (c^ t + (5)/2), (38) 

where c\ and C2 are arbitrary constants. Using the initial conditions a:cm(0) = xq and icm(O) — 0, the solution takes 
the form 

_ xoC(0,a2,(c^t + ^)/2) 
"""-^^^^ - C(0,a2,5/2) ■ ^^^^ 

This solution is plotted in Fig. [TOl where we also plot the result of the numerical solution of Eq. (|36p for the same 
parameters. The agreement between the numerical and the exact result is evident. The advantage of the exact 
analytical solution over the numerical one is that we can investigate the long-time dynamics of the soliton. In Fig. [Til 
we plot the center-of-mass of the soliton for a much longer time interval than in Figs. [S] and [TOl This figure shows 
that the soliton will be trapped over such a large time scale and is oscillating between xo and — a;o- The frequency 
of this oscillation is given by the period of the Mathieu function C(0, a2, (wi -|- (5)/2). A numerical computation of 
the first root of this function for different values of ai and w shows that this frequency is proportional to a^'jJ- The 
constant of proportionality is determined numerically which results in ujp = 0.353 0:2 ^- The fact that the soliton is 
trap ped by the oscillating harmonic potential is a well-established result for such a configuration, known as Paul trap 
[42| . which is used to trap cold ions. (Hence, ujp denotes the frequency of the Paul trap.). 

To have realistic estimates of the parameters ai, a2, and w, we consider the experiment of Strecker et al. [5]. In this 
experiment, solitons were created with a maximum number of iV = 5000 ^Li atoms per soliton. The solitons' center- 
of-mass oscillated with amplitude ~ 370 /xm and period 310 ms. The strength of the harmonic trapping potential in 
the radial direction — 2tt x 800 rad/s was much larger than that of the axial direction ujx — 2tt x 3 rad/s. In this 
case, the unit of length used in this paper is a^; « 2 fim and the unit of time is l/w^, « 50 ms. Furthermore, for the 
^Li scattering length ^ 3ao — 1.5 x 10~^° m, our scaled scattering length is a = as/a^ ~ 10"''. In view of these 
experimental values. Fig. [5] is explained as follows. The filled circles show a soliton located initially at 10 /im from 
the trap center. Oscillating the trapping potential with frequency = 0.44 w^; ~ 27r x 1.3 rad/s, the soliton will be 
drifted a distance of 10 /xm from its initial location in a time period of 1.2 s. On the other hand, using a rather more 
gentle oscillation in the trapping potential, namely with = 0.14 ~ 27r x 0.4 rad/s, the soliton will be nearly 
localized about its initial position over the same time period. 

The exact result, Eq. (|33p . indicates that the soliton can be trapped at any position and with any trap frequency as 
long as the trapping potential and the interatomic interaction are oscillating coherently. Furthermore, it shows that 
the soliton maintains its single-solitonic structure irrespective of the robustness of the trap oscillation. In the present 
case, where the interatomic interactions oscillation is turned off, the situation may be different. The soliton looses 
the possibility of indefinite localization and may also loose its single-soliton structure. Solving Eq. ((36)) for larger 
values of and xq shows that the soliton will leave its initial position faster and gets fragmented into many solitons 
that collide and interfere with each other. Therefore, localizing the soliton at larger distances while maintaining its 
single-solitonic structure requires shallower traps. 
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V. CONCLUSIONS 



Using the Darboux transformation method, we derived exact sohtonc sohitions of a class of Gross-Pitaevskii equa- 
tions represented by Eq. ([7])- The solutions are obtained for a general time-dependent strength of the harmonic 
trapping potential and a related time-dependent strength of the interatomic interaction. Two classes of exact soli- 
tonic solitions were found. The first class represents a single soliton with an arbitrary phase, initial position, and 
initial velocity, as given by Eq. (P^ . The second class comprises single, multiple, and solitary wave solitonic solutions. 

As a specific case, we considered an oscillating trapping potential and interatomic interaction, as given by Eqs. (|30p 
and pip . We found that the soliton can be localized at an arbitrary position and for any amplitude and frequency 
of the oscillating trapping potential. This localization is possible in spite of the asymmetric oscillation where the 
trapping potential spends more time being expulsive. In such a case, one expects that by time the soliton will be 
expelled out of the trap. It turns out, however, that the inhomogeneity in the trapping potential has a balancing 
effect such that it becomes possible to localize the soliton indefinitely. As a consequence of the previously-mentioned 
approximate nature of the Gross-Pitaevskii equation, the fact that our solutions have a ^-function center-of-mass is 
also approximate [43j . For finite number of atoms, the center-of-mass spreads which may lead to delocalization of 
the soliton. For small variations in the center-of-mass, the soliton is expected to remain localized, but when the 
center-of-mass spreading is larger than the amplitude of the soliton oscillation around its equilibrium point, we expect 
that localization disappears completely. 

To discuss the experimental realization, we considered a simpler situation with an oscillating trapping potential and 
constant interatomic interaction, as described by Eq. ([36|) . The numerical solution of this equation showed that soliton 
localization is possible but for a finite time that can be controlled by the frequency and amplitude of the trapping 
potential oscillation. It is shown that for the ^Li experiment of Strecker et al. [5,], the soliton can be localized for a 
time long enough to be observed. With smaller frequency and amplitude, the localization time becomes even larger. 

The Gross-Pitaevskii equation provides an accurate description of the dynamics of solitons as long as finite- 
temperature effects are suppressed and atom losses are negligible (2^ . At finite temperatures and with atom losses, the 
soliton broadens and starts to loose it particle-like behavior. When exposed to oscillations in the trapping potential, 
the soliton will, in this case, be fragmented and soliton localization may not hold. 

Using the localization mechanism, solitons can be prepared in arbitrary initial conditions. This may be useful for 
studying soliton-soliton collisions or the interaction of solitons with potentials. For instance, the shallow oscillating 
trap can be turned on immediately after the solitons were created which leads to localizing the soliton at a certain 
position. Then by switching on and off the oscillating trap the soliton can be moved from one point to the other. 
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APPENDIX A: EXACT SOLUTIONS USING A NONZERO SEED 



A nontrivial seed solution can be easily obtained by substituting in Eq. ([7]) ipix, t) — exp {hi{t) + ih2{x, t)), where 
hi{t) and h2(x,t) are real functions: 



'4'a{x, t) = A exp 



- (4aA2-4e^(*)a;A;o-2/fco' + 4 / e^^^dt {-2 a + ka") + x" ^{t) 



(Al) 



Here A and fco are arbitrary constants. 

Solving the linear system (|4]) and (O using the seed solution, Eq. (|Aip . and then substituting for V'o(a;,t), 'ipi{x,t), 
and 4>i{x,t) in the last equation, we obtain the following new exact solution of Eq. ([7|: 



)(C393 



* 2ari 



2iAy/a: 



ac^e 



/ c^re^^^-^'-v^^ -I- c2re-2°';+A,.v^:'= _ 4AciC2\/^[(2Ai^ - A^)(cos26'i + cos6'3) 

+ (A, - 2Air)(sin26'i +sin6'3) (sin26'i +sin6'3)fco) |, (A2) 
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where 

gi = i[2A'^a{2T] - 1) + ko{ko{l - 2ij) + 2y/7j)x]/2, 

q2 = -2ar] - A^/rjx, 

qs = Ai + ko + i{Ar + 2iXu - 2Xir), 

6*1 = A^Ai^T? - ((2Ah + fco) - Vrj)A,/2, 

92 = (fco(l - 2??) - 4A^Ai,.77)/4 + 2^^0(1 - 2ri) + 2A,(2Ai,?7 " Vvx) + 2fco(A^r] + V^a;), 
6*3 = (2A,Ai, - 2A^Ai,. + A,fco)?7 - ^/riAiX, 
Oi = ko^jx + {2A^a ~ k^){2rj - l)/2, 

r = (A, - 2Ah)2 + (Ar - 2Xir)^ + AA^a + /co(2Ai - 4\u + ko), 
a — A,.Aii + A^Air + Arfco/2, 
A^ = Re[V(2Ai -ifco)2 -4yl2a], 
A, = Im[^(2Ai -ifco)2 -4A2a], 
7 = A, + fco + i{Ar + 2iXi, - 2Xir), 

and ci and C2 are arbitrary constants. It should be noted that this is an exact solution of Eq. ([7]) for any j{t). 

There are 5 arbitrary parameters in the general solution, namely fcp, A, Ai, Ci, and C2. The first three parameters 
control the phase and amplitude of the seed solution which is part of the general solution. The last two parameters 
control the amplitude and phase of the general solution. The solitonic solutions represented by Eq. (|A2p are nonsin- 
gular for all x and t since the denominator of this equation does not vanish. This can be easily deduced from Eq. ([20]) 
where we see that the denominator of Eq. (jA2p is merely the amplitude of 0i and ipi that vanishes only if we have 
the trivial solution with ci = C2 = 0. 
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FIG. 1: Soliton's density profile (solid curve) and trapping potential (dashed curve). First column of subfigures: vo = 
0.8 X xoa2U} sin((5)/2 = 8.4. Middle column of subfigures: vo = xoa2U! sin((5)/2 = 10.5. Last column of subfigures: vo = 
1.2 X xoa2U! sin(5)/2 = 12.6. For the three cases, we take ai = 1, a2 = 0.5, w = 2, 5 = 1, and xo = 50. Time is in units of 
l/wx. 
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FIG. 2: Soliton's center-of-mass trajectory for three values of the initial velocity. The upper curve corresponds to vo = 
1.2 X xo a2U) sin {5)/2 = 12.6, the middle curve corresponds to vo = xo aiuj sin ((5)/2 = 10.5, and the lower curve corresponds 
to Vo = 0.8 X a;o 0:20; sin ((5)/2 = 8.4. For the three curves, we take ai = \, ai = 0.5, uj = 2, 6 = \, and xo = 50. 




FIG. 3: Strength of the trapping potential versus time. The following values were used: a\ = 1, a2 = 0.75, a; = tt, 5 = 1.2. 
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FIG. 4: First column: The trapping potential V{xcin{t)). Second column: The effective trapping potential V{xcin{t)). (a) 
and (b) correspond to, vo = 1.2a2XQ0j sin((5)/2, (c) and (d) correspond to, vo =02X00; sin(5)/2, (e) and (f) correspond to, 
vo = 0.80:2 xq uj sin {5)/2. For all curves: xo = 50, oi = 1, 02 = 0.2, 5 = 0.4, and w = tt. The arrows show the direction of the 
motion of the soliton. 



FIG. 5: (Left) Spatiotemporal contour plot of soliton density profile. (Right) Soliton density profile at t = 0. The values of 
the parameters used in this plot are: ci = — C2 = 10, 03=04 = 1, Aij = Air = 1, A = 2, a = 0.9, fco = 5, A = a; = 1, 5 = 0, 
ai = —6, a2 = 0.3. The value of ko is the solution of a = with respect to fco. 
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FIG. 7: (Color online) Density profile of a single soliton solution. The values of the parameters used in this plot are: — C2 = 
ci — 10, C3 = C4 = 1, Aii — Air — 1, A = 2, a = 0.9, ko ~ 5.03, X = lo = 1, 6 = 0, ai = —6, 02 = 0.3. The value of ko is the 
solution of a = with respect to ko. 
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FIG. 8: (Color online) Density profile of a multi-solitons solution. The values of the parameters used in this plot are: — C2 = 
ci = 10, C3 = C4 = 1, Aii = 0, Ai^ = 1, ^ = 2, a = 0.9, fco = 0, A = 1, a; = 0.01, 5 = 0, ai = -2, 02 = 0.3. 
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FIG. 9: Soliton's center-of-mass trajectory. The solid and dashed curves correspond to the exact formula Eq. (I24|l . Empty and 
filled circles correspond to the numerical solution of Eq. (|36p . Solid curve and filled circles are obtained with Q = LOyUi2/2 = 
0.44: ujx and dashed curve and empty circles are obtained with Q — oj-^/aT/S = 0.14 uj^- The rest of parameters used are: 
a = 10"*, 5 = ?;o = 0, xo = 5, ai = 1, A'^ = 4 X 10^, 02 — 0.1 for the solid curve and 02 = 0.01 for the dashed curve. 
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FIG. 11: The solid and dashed curves correspond to the exact formula Eq. (|39p . The two curves are the same as those of Fig. 1101 
but shown here over a larger time interval. Extending the time interval further shows that the dashed curve is also oscillatory. 
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